############################################################################
rm(list=ls())
library(readxl)
library(readr)
library(countrycode)
library(QuantPsyc)
###################################################################################################
#            Regression Plots for Constitutional Rigidity and Amendment Frequency                 #
###################################################################################################
library(ggplot2)
library(readr)

mydata <- read.csv("amendment_table_regression_data_updated.csv")

# 1. All Amendments

myDF <- mydata
myDF$x <- mydata$VP_rigid_eu
myDF$y <- mydata$all_am_rate

#mod_lin <- lm(y~x, data=myDF)
#w <- 1/fitted( lm(abs(residuals(mod_lin))~fitted(mod_lin)) )^2
#t <- lm(y~x, data=myDF, weights=w)

w = 1/myDF$x**2
t = lm(y~x, data=myDF, weights=1/w) 
summary(t)

t.standardized <- lm(scale(y)~scale(x), data=myDF, weights=1/w)
summary(t.standardized)

newdata = data.frame(x=seq(0.5,1.75,0.01))
w = 1/newdata$x**2
p1 = predict.lm(t, newdata, interval = 'prediction', weights=1/w)
a <- ggplot()
a <- a + geom_point(data=myDF, aes(x=x, y=y), shape=1)+geom_text(data=myDF, aes(x=x, y=y,label=abbrev),hjust=0, vjust=0,size=6)
a <- a + geom_abline(intercept=t$coefficients[1],   slope=t$coefficients[2], color='blue')  
newdata$lwr = p1[,c("lwr")]
newdata$upr = p1[,c("upr")]
a <- a #+ geom_ribbon(data=newdata, aes(x=x, ymin=lwr, ymax=upr),   fill='gray', alpha=0.5)
a <- a + ggtitle("VP Rigidity and Amendment Frequency") + xlab("Veto Player Rigidity") + ylab("Frequency of All Amendments")
a <- a + theme(
  plot.title = element_text(color = "black", size = 20, face = "bold"),
  axis.title.x = element_text(size = 15),
  axis.title.y = element_text(size = 15)
)
a <- a+annotate("text", x = 0.55, y = 1, label = "Coefficient= -0.33636") +
      annotate("text", x = 0.55, y = 0.975, label = "P-Value = 0.00032") 

a

